home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / odepack / stode.f < prev    next >
Text File  |  1995-12-14  |  20KB  |  484 lines
  1.       SUBROUTINE STODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR,
  2.      1   WM, IWM, F, JAC, PJAC, SLVS, IERR)
  3. CLLL. OPTIMIZE
  4.       EXTERNAL F, JAC, PJAC, SLVS
  5.       INTEGER NEQ, NYH, IWM
  6.       INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
  7.      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
  8.      2   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  9.       INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ
  10.       DOUBLE PRECISION Y, YH, YH1, EWT, SAVF, ACOR, WM
  11.       DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
  12.      2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
  13.       DOUBLE PRECISION DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP,
  14.      1   R, RH, RHDN, RHSM, RHUP, TOLD, VNORM
  15.       DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*),
  16.      1   ACOR(*), WM(*), IWM(*)
  17.       COMMON /LS0001/ CONIT, CRATE, EL(13), ELCO(13,12),
  18.      1   HOLD, RMAX, TESCO(3,12),
  19.      2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, IOWND(14),
  20.      3   IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
  21.      4   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
  22.      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  23. C-----------------------------------------------------------------------
  24. C STODE PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE
  25. C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS.
  26. C NOTE.. STODE IS INDEPENDENT OF THE VALUE OF THE ITERATION METHOD
  27. C INDICATOR MITER, WHEN THIS IS .NE. 0, AND HENCE IS INDEPENDENT
  28. C OF THE TYPE OF CHORD METHOD USED, OR THE JACOBIAN STRUCTURE.
  29. C COMMUNICATION WITH STODE IS DONE WITH THE FOLLOWING VARIABLES..
  30. C
  31. C NEQ    = INTEGER ARRAY CONTAINING PROBLEM SIZE IN NEQ(1), AND
  32. C          PASSED AS THE NEQ ARGUMENT IN ALL CALLS TO F AND JAC.
  33. C Y      = AN ARRAY OF LENGTH .GE. N USED AS THE Y ARGUMENT IN
  34. C          ALL CALLS TO F AND JAC.
  35. C YH     = AN NYH BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES
  36. C          AND THEIR APPROXIMATE SCALED DERIVATIVES, WHERE
  37. C          LMAX = MAXORD + 1.  YH(I,J+1) CONTAINS THE APPROXIMATE
  38. C          J-TH DERIVATIVE OF Y(I), SCALED BY H**J/FACTORIAL(J)
  39. C          (J = 0,1,...,NQ).  ON ENTRY FOR THE FIRST STEP, THE FIRST
  40. C          TWO COLUMNS OF YH MUST BE SET FROM THE INITIAL VALUES.
  41. C NYH    = A CONSTANT INTEGER .GE. N, THE FIRST DIMENSION OF YH.
  42. C YH1    = A ONE-DIMENSIONAL ARRAY OCCUPYING THE SAME SPACE AS YH.
  43. C EWT    = AN ARRAY OF LENGTH N CONTAINING MULTIPLICATIVE WEIGHTS
  44. C          FOR LOCAL ERROR MEASUREMENTS.  LOCAL ERRORS IN Y(I) ARE
  45. C          COMPARED TO 1.0/EWT(I) IN VARIOUS ERROR TESTS.
  46. C SAVF   = AN ARRAY OF WORKING STORAGE, OF LENGTH N.
  47. C          ALSO USED FOR INPUT OF YH(*,MAXORD+2) WHEN JSTART = -1
  48. C          AND MAXORD .LT. THE CURRENT ORDER NQ.
  49. C ACOR   = A WORK ARRAY OF LENGTH N, USED FOR THE ACCUMULATED
  50. C          CORRECTIONS.  ON A SUCCESSFUL RETURN, ACOR(I) CONTAINS
  51. C          THE ESTIMATED ONE-STEP LOCAL ERROR IN Y(I).
  52. C WM,IWM = REAL AND INTEGER WORK ARRAYS ASSOCIATED WITH MATRIX
  53. C          OPERATIONS IN CHORD ITERATION (MITER .NE. 0).
  54. C PJAC   = NAME OF ROUTINE TO EVALUATE AND PREPROCESS JACOBIAN MATRIX 
  55. C          AND P = I - H*EL0*JAC, IF A CHORD METHOD IS BEING USED.
  56. C SLVS   = NAME OF ROUTINE TO SOLVE LINEAR SYSTEM IN CHORD ITERATION. 
  57. C CCMAX  = MAXIMUM RELATIVE CHANGE IN H*EL0 BEFORE PJAC IS CALLED.
  58. C H      = THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP.
  59. C          H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE
  60. C          PROBLEM.  H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS
  61. C          SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM.
  62. C HMIN   = THE MINIMUM ABSOLUTE VALUE OF THE STEP SIZE H TO BE USED.
  63. C HMXI   = INVERSE OF THE MAXIMUM ABSOLUTE VALUE OF H TO BE USED.
  64. C          HMXI = 0.0 IS ALLOWED AND CORRESPONDS TO AN INFINITE HMAX. 
  65. C          HMIN AND HMXI MAY BE CHANGED AT ANY TIME, BUT WILL NOT
  66. C          TAKE EFFECT UNTIL THE NEXT CHANGE OF H IS CONSIDERED.
  67. C TN     = THE INDEPENDENT VARIABLE. TN IS UPDATED ON EACH STEP TAKEN.
  68. C JSTART = AN INTEGER USED FOR INPUT ONLY, WITH THE FOLLOWING
  69. C          VALUES AND MEANINGS..
  70. C               0  PERFORM THE FIRST STEP.
  71. C           .GT.0  TAKE A NEW STEP CONTINUING FROM THE LAST.
  72. C              -1  TAKE THE NEXT STEP WITH A NEW VALUE OF H, MAXORD,
  73. C                    N, METH, MITER, AND/OR MATRIX PARAMETERS.
  74. C              -2  TAKE THE NEXT STEP WITH A NEW VALUE OF H,
  75. C                    BUT WITH OTHER INPUTS UNCHANGED.
  76. C          ON RETURN, JSTART IS SET TO 1 TO FACILITATE CONTINUATION.
  77. C KFLAG  = A COMPLETION CODE WITH THE FOLLOWING MEANINGS..
  78. C               0  THE STEP WAS SUCCESFUL.
  79. C              -1  THE REQUESTED ERROR COULD NOT BE ACHIEVED.
  80. C              -2  CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED.
  81. C              -3  FATAL ERROR IN PJAC OR SLVS.
  82. C          A RETURN WITH KFLAG = -1 OR -2 MEANS EITHER
  83. C          ABS(H) = HMIN OR 10 CONSECUTIVE FAILURES OCCURRED.
  84. C          ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF TN AND
  85. C          THE YH ARRAY ARE AS OF THE BEGINNING OF THE LAST 
  86. C          STEP, AND H IS THE LAST STEP SIZE ATTEMPTED.
  87. C MAXORD = THE MAXIMUM ORDER OF INTEGRATION METHOD TO BE ALLOWED.
  88. C MAXCOR = THE MAXIMUM NUMBER OF CORRECTOR ITERATIONS ALLOWED.
  89. C MSBP   = MAXIMUM NUMBER OF STEPS BETWEEN PJAC CALLS (MITER .GT. 0). 
  90. C MXNCF  = MAXIMUM NUMBER OF CONVERGENCE FAILURES ALLOWED.
  91. C METH/MITER = THE METHOD FLAGS.  SEE DESCRIPTION IN DRIVER.
  92. C N      = THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS.
  93. C IERR   = ERROR FLAG FROM USER-SUPPLIED FUNCTION
  94. C-----------------------------------------------------------------------
  95.       KFLAG = 0
  96.       TOLD = TN
  97.       NCF = 0
  98.       IERPJ = 0
  99.       IERSL = 0
  100.       JCUR = 0
  101.       ICF = 0
  102.       DELP = 0.0D0
  103.       IF (JSTART .GT. 0) GO TO 200
  104.       IF (JSTART .EQ. -1) GO TO 100
  105.       IF (JSTART .EQ. -2) GO TO 160
  106. C-----------------------------------------------------------------------
  107. C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER VARIABLES ARE
  108. C INITIALIZED.  RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED 
  109. C IN A SINGLE STEP.  IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL 
  110. C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10.  IF A FAILURE
  111. C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2
  112. C FOR THE NEXT INCREASE.
  113. C-----------------------------------------------------------------------
  114.       LMAX = MAXORD + 1
  115.       NQ = 1
  116.       L = 2
  117.       IALTH = 2
  118.       RMAX = 10000.0D0
  119.       RC = 0.0D0
  120.       EL0 = 1.0D0
  121.       CRATE = 0.7D0 
  122.       HOLD = H
  123.       MEO = METH
  124.       NSLP = 0
  125.       IPUP = MITER
  126.       IRET = 3
  127.       GO TO 140
  128. C-----------------------------------------------------------------------
  129. C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN JSTART = -1.
  130. C IPUP IS SET TO MITER TO FORCE A MATRIX UPDATE.
  131. C IF AN ORDER INCREASE IS ABOUT TO BE CONSIDERED (IALTH = 1),
  132. C IALTH IS RESET TO 2 TO POSTPONE CONSIDERATION ONE MORE STEP.
  133. C IF THE CALLER HAS CHANGED METH, CFODE IS CALLED TO RESET
  134. C THE COEFFICIENTS OF THE METHOD.
  135. C IF THE CALLER HAS CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT
  136. C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN ACCORDINGLY.
  137. C IF H IS TO BE CHANGED, YH MUST BE RESCALED.
  138. C IF H OR METH IS BEING CHANGED, IALTH IS RESET TO L = NQ + 1
  139. C TO PREVENT FURTHER CHANGES IN H FOR THAT MANY STEPS.
  140. C-----------------------------------------------------------------------
  141.  100  IPUP = MITER
  142.       LMAX = MAXORD + 1
  143.       IF (IALTH .EQ. 1) IALTH = 2
  144.       IF (METH .EQ. MEO) GO TO 110
  145.       CALL CFODE (METH, ELCO, TESCO)
  146.       MEO = METH
  147.       IF (NQ .GT. MAXORD) GO TO 120
  148.       IALTH = L
  149.       IRET = 1
  150.       GO TO 150
  151.  110  IF (NQ .LE. MAXORD) GO TO 160
  152.  120  NQ = MAXORD
  153.       L = LMAX
  154.       DO 125 I = 1,L
  155.  125    EL(I) = ELCO(I,NQ)
  156.       NQNYH = NQ*NYH
  157.       RC = RC*EL(1)/EL0
  158.       EL0 = EL(1)
  159.       CONIT = 0.5D0/DBLE(NQ+2)
  160.       DDN = VNORM (N, SAVF, EWT)/TESCO(1,L)
  161.       EXDN = 1.0D0/DBLE(L)  
  162.       RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
  163.       RH = DMIN1(RHDN,1.0D0)
  164.       IREDO = 3
  165.       IF (H .EQ. HOLD) GO TO 170
  166.       RH = DMIN1(RH,DABS(H/HOLD))
  167.       H = HOLD
  168.       GO TO 175
  169. C-----------------------------------------------------------------------
  170. C CFODE IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS FOR THE
  171. C CURRENT METH.  THEN THE EL VECTOR AND RELATED CONSTANTS ARE RESET
  172. C WHENEVER THE ORDER NQ IS CHANGED, OR AT THE START OF THE PROBLEM.
  173. C-----------------------------------------------------------------------
  174.  140  CALL CFODE (METH, ELCO, TESCO)
  175.  150  DO 155 I = 1,L
  176.  155    EL(I) = ELCO(I,NQ)
  177.       NQNYH = NQ*NYH
  178.       RC = RC*EL(1)/EL0
  179.       EL0 = EL(1)
  180.       CONIT = 0.5D0/DBLE(NQ+2)
  181.       GO TO (160, 170, 200), IRET
  182. C-----------------------------------------------------------------------
  183. C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST
  184. C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED.  IALTH IS SET TO
  185. C L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT MANY STEPS, UNLESS
  186. C FORCED BY A CONVERGENCE OR ERROR TEST FAILURE.
  187. C-----------------------------------------------------------------------
  188.  160  IF (H .EQ. HOLD) GO TO 200
  189.       RH = H/HOLD
  190.       H = HOLD
  191.       IREDO = 3
  192.       GO TO 175
  193.  170  RH = DMAX1(RH,HMIN/DABS(H))
  194.  175  RH = DMIN1(RH,RMAX)
  195.       RH = RH/DMAX1(1.0D0,DABS(H)*HMXI*RH)
  196.       R = 1.0D0
  197.       DO 180 J = 2,L
  198.         R = R*RH
  199.         DO 180 I = 1,N
  200.  180      YH(I,J) = YH(I,J)*R 
  201.       H = H*RH
  202.       RC = RC*RH
  203.       IALTH = L
  204.       IF (IREDO .EQ. 0) GO TO 690
  205. C-----------------------------------------------------------------------
  206. C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY 
  207. C MULTIPLYING THE YH ARRAY BY THE PASCAL TRIANGLE MATRIX.
  208. C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT  H*EL(1).
  209. C WHEN RC DIFFERS FROM 1 BY MORE THAN CCMAX, IPUP IS SET TO MITER
  210. C TO FORCE PJAC TO BE CALLED, IF A JACOBIAN IS INVOLVED.
  211. C IN ANY CASE, PJAC IS CALLED AT LEAST EVERY MSBP STEPS.
  212. C-----------------------------------------------------------------------
  213.  200  IF (DABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER 
  214.       IF (NST .GE. NSLP+MSBP) IPUP = MITER
  215.       TN = TN + H
  216.       I1 = NQNYH + 1
  217.       DO 215 JB = 1,NQ
  218.         I1 = I1 - NYH
  219. CDIR$ IVDEP
  220.         DO 210 I = I1,NQNYH
  221.  210      YH1(I) = YH1(I) + YH1(I+NYH)
  222.  215    CONTINUE
  223. C-----------------------------------------------------------------------
  224. C UP TO MAXCOR CORRECTOR ITERATIONS ARE TAKEN.  A CONVERGENCE TEST IS 
  225. C MADE ON THE R.M.S. NORM OF EACH CORRECTION, WEIGHTED BY THE ERROR
  226. C WEIGHT VECTOR EWT.  THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE
  227. C VECTOR ACOR(I).  THE YH ARRAY IS NOT ALTERED IN THE CORRECTOR LOOP. 
  228. C-----------------------------------------------------------------------
  229.  220  M = 0
  230.       DO 230 I = 1,N
  231.  230    Y(I) = YH(I,1)
  232.       IERR = 0
  233.       CALL F (NEQ, TN, Y, SAVF, IERR)
  234.       IF (IERR .LT. 0) RETURN
  235.       NFE = NFE + 1 
  236.       IF (IPUP .LE. 0) GO TO 250
  237. C-----------------------------------------------------------------------
  238. C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED AND
  239. C PREPROCESSED BEFORE STARTING THE CORRECTOR ITERATION.  IPUP IS SET
  240. C TO 0 AS AN INDICATOR THAT THIS HAS BEEN DONE.
  241. C-----------------------------------------------------------------------
  242.       IERR = 0
  243.       CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC,
  244.      1   IERR)
  245.       IF (IERR .LT. 0) RETURN
  246.       IPUP = 0
  247.       RC = 1.0D0
  248.       NSLP = NST
  249.       CRATE = 0.7D0 
  250.       IF (IERPJ .NE. 0) GO TO 430
  251.  250  DO 260 I = 1,N
  252.  260    ACOR(I) = 0.0D0
  253.  270  IF (MITER .NE. 0) GO TO 350
  254. C-----------------------------------------------------------------------
  255. C IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM
  256. C THE RESULT OF THE LAST FUNCTION EVALUATION.
  257. C-----------------------------------------------------------------------
  258.       DO 290 I = 1,N
  259.         SAVF(I) = H*SAVF(I) - YH(I,2)
  260.  290    Y(I) = SAVF(I) - ACOR(I)
  261.       DEL = VNORM (N, Y, EWT) 
  262.       DO 300 I = 1,N
  263.         Y(I) = YH(I,1) + EL(1)*SAVF(I)
  264.  300    ACOR(I) = SAVF(I)
  265.       GO TO 400
  266. C-----------------------------------------------------------------------
  267. C IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR,
  268. C AND SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND SIDE AND
  269. C P AS COEFFICIENT MATRIX.
  270. C-----------------------------------------------------------------------
  271.  350  DO 360 I = 1,N
  272.  360    Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
  273.       CALL SLVS (WM, IWM, Y, SAVF)
  274.       IF (IERSL .LT. 0) GO TO 430
  275.       IF (IERSL .GT. 0) GO TO 410
  276.       DEL = VNORM (N, Y, EWT) 
  277.       DO 380 I = 1,N
  278.         ACOR(I) = ACOR(I) + Y(I)
  279.  380    Y(I) = YH(I,1) + EL(1)*ACOR(I)
  280. C-----------------------------------------------------------------------
  281. C TEST FOR CONVERGENCE.  IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE
  282. C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST.
  283. C-----------------------------------------------------------------------
  284.  400  IF (M .NE. 0) CRATE = DMAX1(0.2D0*CRATE,DEL/DELP)
  285.       DCON = DEL*DMIN1(1.0D0,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT)
  286.       IF (DCON .LE. 1.0D0) GO TO 450
  287.       M = M + 1
  288.       IF (M .EQ. MAXCOR) GO TO 410
  289.       IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410
  290.       DELP = DEL
  291.       IERR = 0
  292.       CALL F (NEQ, TN, Y, SAVF, IERR)
  293.       IF (IERR .LT. 0) RETURN
  294.       NFE = NFE + 1 
  295.       GO TO 270
  296. C-----------------------------------------------------------------------
  297. C THE CORRECTOR ITERATION FAILED TO CONVERGE.
  298. C IF MITER .NE. 0 AND THE JACOBIAN IS OUT OF DATE, PJAC IS CALLED FOR 
  299. C THE NEXT TRY.  OTHERWISE THE YH ARRAY IS RETRACTED TO ITS VALUES
  300. C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE.  IF H CANNOT BE
  301. C REDUCED OR MXNCF FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -2.
  302. C-----------------------------------------------------------------------
  303.  410  IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430
  304.       ICF = 1
  305.       IPUP = MITER
  306.       GO TO 220
  307.  430  ICF = 2
  308.       NCF = NCF + 1 
  309.       RMAX = 2.0D0
  310.       TN = TOLD
  311.       I1 = NQNYH + 1
  312.       DO 445 JB = 1,NQ
  313.         I1 = I1 - NYH
  314. CDIR$ IVDEP
  315.         DO 440 I = I1,NQNYH
  316.  440      YH1(I) = YH1(I) - YH1(I+NYH)
  317.  445    CONTINUE
  318.       IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680
  319.       IF (DABS(H) .LE. HMIN*1.00001D0) GO TO 670
  320.       IF (NCF .EQ. MXNCF) GO TO 670
  321.       RH = 0.25D0
  322.       IPUP = MITER
  323.       IREDO = 1
  324.       GO TO 170
  325. C-----------------------------------------------------------------------
  326. C THE CORRECTOR HAS CONVERGED.  JCUR IS SET TO 0
  327. C TO SIGNAL THAT THE JACOBIAN INVOLVED MAY NEED UPDATING LATER.
  328. C THE LOCAL ERROR TEST IS MADE AND CONTROL PASSES TO STATEMENT 500
  329. C IF IT FAILS.
  330. C-----------------------------------------------------------------------
  331.  450  JCUR = 0
  332.       IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
  333.       IF (M .GT. 0) DSM = VNORM (N, ACOR, EWT)/TESCO(2,NQ)
  334.       IF (DSM .GT. 1.0D0) GO TO 500
  335. C-----------------------------------------------------------------------
  336. C AFTER A SUCCESSFUL STEP, UPDATE THE YH ARRAY.
  337. C CONSIDER CHANGING H IF IALTH = 1.  OTHERWISE DECREASE IALTH BY 1.
  338. C IF IALTH IS THEN 1 AND NQ .LT. MAXORD, THEN ACOR IS SAVED FOR
  339. C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP.
  340. C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER
  341. C BY ONE IS CONSIDERED ALSO.  A CHANGE IN H IS MADE ONLY IF IT IS BY A
  342. C FACTOR OF AT LEAST 1.1.  IF NOT, IALTH IS SET TO 3 TO PREVENT
  343. C TESTING FOR THAT MANY STEPS.
  344. C-----------------------------------------------------------------------
  345.       KFLAG = 0
  346.       IREDO = 0
  347.       NST = NST + 1 
  348.       HU = H
  349.       NQU = NQ
  350.       DO 470 J = 1,L
  351.         DO 470 I = 1,N
  352.  470      YH(I,J) = YH(I,J) + EL(J)*ACOR(I)
  353.       IALTH = IALTH - 1
  354.       IF (IALTH .EQ. 0) GO TO 520
  355.       IF (IALTH .GT. 1) GO TO 700
  356.       IF (L .EQ. LMAX) GO TO 700
  357.       DO 490 I = 1,N
  358.  490    YH(I,LMAX) = ACOR(I)
  359.       GO TO 700
  360. C-----------------------------------------------------------------------
  361. C THE ERROR TEST FAILED.  KFLAG KEEPS TRACK OF MULTIPLE FAILURES.
  362. C RESTORE TN AND THE YH ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE
  363. C TO TRY THE STEP AGAIN.  COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR
  364. C ONE LOWER ORDER.  AFTER 2 OR MORE FAILURES, H IS FORCED TO DECREASE 
  365. C BY A FACTOR OF 0.2 OR LESS. 
  366. C-----------------------------------------------------------------------
  367.  500  KFLAG = KFLAG - 1
  368.       TN = TOLD
  369.       I1 = NQNYH + 1
  370.       DO 515 JB = 1,NQ
  371.         I1 = I1 - NYH
  372. CDIR$ IVDEP
  373.         DO 510 I = I1,NQNYH
  374.  510      YH1(I) = YH1(I) - YH1(I+NYH)
  375.  515    CONTINUE
  376.       RMAX = 2.0D0
  377.       IF (DABS(H) .LE. HMIN*1.00001D0) GO TO 660
  378.       IF (KFLAG .LE. -3) GO TO 640
  379.       IREDO = 2
  380.       RHUP = 0.0D0
  381.       GO TO 540
  382. C-----------------------------------------------------------------------
  383. C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS 
  384. C RHDN, RHSM, AND RHUP ARE COMPUTED, BY WHICH H COULD BE MULTIPLIED
  385. C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY. 
  386. C IN THE CASE OF FAILURE, RHUP = 0.0 TO AVOID AN ORDER INCREASE.
  387. C THE LARGEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN
  388. C ACCORDINGLY.  IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE
  389. C ADDITIONAL SCALED DERIVATIVE.
  390. C-----------------------------------------------------------------------
  391.  520  RHUP = 0.0D0
  392.       IF (L .EQ. LMAX) GO TO 540
  393.       DO 530 I = 1,N
  394.  530    SAVF(I) = ACOR(I) - YH(I,LMAX)
  395.       DUP = VNORM (N, SAVF, EWT)/TESCO(3,NQ)
  396.       EXUP = 1.0D0/DBLE(L+1)
  397.       RHUP = 1.0D0/(1.4D0*DUP**EXUP + 0.0000014D0)
  398.  540  EXSM = 1.0D0/DBLE(L)  
  399.       RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0)
  400.       RHDN = 0.0D0
  401.       IF (NQ .EQ. 1) GO TO 560
  402.       DDN = VNORM (N, YH(1,L), EWT)/TESCO(1,NQ)
  403.       EXDN = 1.0D0/DBLE(NQ) 
  404.       RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
  405.  560  IF (RHSM .GE. RHUP) GO TO 570
  406.       IF (RHUP .GT. RHDN) GO TO 590
  407.       GO TO 580
  408.  570  IF (RHSM .LT. RHDN) GO TO 580
  409.       NEWQ = NQ
  410.       RH = RHSM
  411.       GO TO 620
  412.  580  NEWQ = NQ - 1 
  413.       RH = RHDN
  414.       IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0
  415.       GO TO 620
  416.  590  NEWQ = L
  417.       RH = RHUP
  418.       IF (RH .LT. 1.1D0) GO TO 610
  419.       R = EL(L)/DBLE(L)     
  420.       DO 600 I = 1,N
  421.  600    YH(I,NEWQ+1) = ACOR(I)*R
  422.       GO TO 630
  423.  610  IALTH = 3
  424.       GO TO 700
  425.  620  IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1D0)) GO TO 610
  426.       IF (KFLAG .LE. -2) RH = DMIN1(RH,0.2D0)
  427. C-----------------------------------------------------------------------
  428. C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS.
  429. C IN ANY CASE H IS RESET ACCORDING TO RH AND THE YH ARRAY IS RESCALED.
  430. C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE.
  431. C-----------------------------------------------------------------------
  432.       IF (NEWQ .EQ. NQ) GO TO 170
  433.  630  NQ = NEWQ
  434.       L = NQ + 1
  435.       IRET = 2
  436.       GO TO 150
  437. C-----------------------------------------------------------------------
  438. C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURED.
  439. C IF 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -1.
  440. C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE
  441. C YH ARRAY HAVE ERRORS OF THE WRONG ORDER.  HENCE THE FIRST 
  442. C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1.  THEN
  443. C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED,
  444. C UNTIL IT SUCCEEDS OR H REACHES HMIN.
  445. C-----------------------------------------------------------------------
  446.  640  IF (KFLAG .EQ. -10) GO TO 660
  447.       RH = 0.1D0
  448.       RH = DMAX1(HMIN/DABS(H),RH)
  449.       H = H*RH
  450.       DO 645 I = 1,N
  451.  645    Y(I) = YH(I,1)
  452.       IERR = 0
  453.       CALL F (NEQ, TN, Y, SAVF, IERR)
  454.       IF (IERR .LT. 0) RETURN
  455.       NFE = NFE + 1 
  456.       DO 650 I = 1,N
  457.  650    YH(I,2) = H*SAVF(I)
  458.       IPUP = MITER
  459.       IALTH = 5
  460.       IF (NQ .EQ. 1) GO TO 200
  461.       NQ = 1
  462.       L = 2
  463.       IRET = 3
  464.       GO TO 150
  465. C-----------------------------------------------------------------------
  466. C ALL RETURNS ARE MADE THROUGH THIS SECTION.  H IS SAVED IN HOLD
  467. C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP.
  468. C-----------------------------------------------------------------------
  469.  660  KFLAG = -1
  470.       GO TO 720
  471.  670  KFLAG = -2
  472.       GO TO 720
  473.  680  KFLAG = -3
  474.       GO TO 720
  475.  690  RMAX = 10.0D0 
  476.  700  R = 1.0D0/TESCO(2,NQU)
  477.       DO 710 I = 1,N
  478.  710    ACOR(I) = ACOR(I)*R
  479.  720  HOLD = H
  480.       JSTART = 1
  481.       RETURN
  482. C----------------------- END OF SUBROUTINE STODE -----------------------
  483.       END 
  484.